#include "solver_sav_onp.h"

#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <stdio.h>
#include <omp.h>
#include "portable_time.h"

double f_test(double num1, double num2)
{
	return rand() / 100;
}

void waitAllThreads(int threadID,int maxThreads, bool threadReady)
{
	do
	{
		bool exit = true;
		for (int i = 0; i < maxThreads; i++)
		{
			if (threadReady == false) exit = false;
		}

	}while (exit == false);
}

double* solver_omp(double* u0, Area area, double epsilon, TwoDimFunc f)
{
	int n = area.n;
	double a = area.a, b = area.b;
	double* u_cur = new double[n * n];
	double* u_prev = new double[n * n];
	double h = (b - a) / n;
	double error = 0;
	double start, finish;

	double* f_cur = new double[n * n];
	int maxThreads = omp_get_max_threads();
	printf("\nThreads = %d\n", maxThreads);
	double* threadError = new double[maxThreads];
	bool* threadReady = new bool[maxThreads];
	start = PortableGetTime();
	omp_set_num_threads(maxThreads);
#pragma omp parallel for
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < n; ++j)
			{
				//int threadID = omp_get_thread_num();
				f_cur[i * n + j] = f(a + i * h, a + j * h);
				//printf("%d %d form Thread %d\n",i,j,threadID);
			}
		}
	finish = PortableGetTime();
	printf("Compute function time: %.4f\n", finish - start);

	memcpy(u_prev, u0, n * n * sizeof(double));
	memcpy(u_cur, u0, n * n * sizeof(double));
	printf("\nu_prev and u_cur created\n\n");
	start = PortableGetTime();
	omp_set_num_threads(maxThreads);
#pragma omp parallel 
{
	int threadID = omp_get_thread_num();
	do
	{ 
		for (int i = threadID + 1; i < n - 1; i+=maxThreads)
		{
			for (int j = 1; j < n - 1; ++j)
			{
				u_cur[i * n + j] = 0.25 * 
					(u_prev[(i - 1) * n + j] + u_prev[(i + 1) * n + j] +
					u_prev[i * n + j - 1] + u_prev[i * n + j + 1] - 
					h * h * f_cur[i * n + j]);
			}
		}
		threadReady[threadID] =	true;
		waitAllThreads(threadID,maxThreads,threadReady);
		threadReady[threadID] = false;

		threadError[threadID] = abs(u_cur[(threadID + 1) * n + 1] - u_prev[(threadID + 1) * n + 1]);
		double e;
		for (int i = threadID + 1; i < n - 1; i+=maxThreads)
		{
			for (int j = 1; j < n - 1; ++j)
			{
				e = abs(u_cur[i*n +j] - u_prev[i*n +j]);
				if (e > threadError[threadID]) threadError[threadID] = e;
			}
		}

		for (int i=0;i<maxThreads;i++)
		{
			if (error < threadError[i]) error = threadError[i];
		}
	
		threadReady[threadID] = true;
		waitAllThreads(threadID,maxThreads,threadReady);
		threadReady[threadID] = false;


	if (threadID == 0)
	{
		error = 0;
		for (int i=0;i<maxThreads;i++)
		{
			if (error < threadError[i]) error = threadError[i];
		}
		printf("error = %f\n", error);
		double* tmp;
		tmp = u_cur; u_cur = u_prev; u_prev = tmp;
	}
	threadReady[threadID] = true;
	waitAllThreads(threadID,maxThreads,threadReady);
	threadReady[threadID] = false;

	} while (error > epsilon);
}
	finish = PortableGetTime();
	printf("Solve time: %.4f\n", finish - start);

	delete[] u_cur;
	return u_prev;
}